Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 23 August 2011 (MN IMfeX style file v2.2) 



Time-Dependent Behavior of Lyman a Photon Transfer in 
High Redshift Optically Thick Medium 



O ■ Wen Xu 1 ' 2 , Xiang-Ping Wu 1 , and Li-Zhi Fang 2 

£\J i 1 National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012, China 
' 2 Department of Physics, University of Arizona, Tucson, AZ 85721, USA 



23 August 2011 



ABSTRACT 

With Monte Carlo simulation method, we investigate the time dependent behavior 
of Lya photon transfer in optically thick medium of the concordance ACDM universe. 
At high redshift, the Lya photon escaping from optically thick medium has a time 
scale as long as the age of the luminous object, or even comparable to the age of the 
universe. In this case, time-independent, or stationary solutions of the Lya photon 
transfer with resonant scattering will overlook important features of the escaped Lya 
photons in physical and frequency spaces. More seriously, the expansion of the universe 
leads to that the time- independent solutions of the Lya photon transfer may not exist. 
We show that time-dependent solutions sometimes are essential for understanding the 
Lya emission and absorption at high redshifts. For Lya photons from sources at 
redshift 1 + z = 10 and being surrounded by neutral hydrogen IGM of the ACDM 
universe, the escape coefficient is found to be always less, or much less than one, 
regardless of the age or life time of the sources. Under such environment, we also find 
that even when the Lya photon luminosity of the sources is stable, the mean surface 
brightness is gradually increasing in the first 10 6 years, and then decreasing with a 
power law of time, but never approaches a stable, time-independent state. That is, 
all 1 + z = 10 sources in a neutral Hubble expanding IGM with Lya luminosity L 
have their maximum of mean surface brightness ~ 10~ 21 (L/10 43 erg/s) erg s _1 cm~ 2 
arcsec -2 at the age of about 10 6 years. The time-dependent effects on the red damping 
wing profile are also addressed. 
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1 INTRODUCTION 



Lya photons have been widely applied to study the physics 
of the universe at redshifts from 2 to 8. The redshifted Lya 
photons carry the information of the photon source, the halo 
surrounding the source, and the IGM at early universe. The 
optical afterglow of gamma-ray burst (GRB) has been mod- 
eled as the red wing of Lya photon absorption, and used 
to estimate the column number density of neutral hydrogen 
of IG M at high redshift (|Totani et al.ll2006l ; ISalvaterra et ail 
2009). The transmitted flux of QSO absorption spectrum at 
redshift z > 5 consists of complete absorption troughs sep - 
arated by spikes (e.g.. iBecker et al.|[200lf ; iFan et aL| [2006). 



The spikes have been explained as Lya photons leaki ng at 
low density areas fe.g. lLiu et aT]|2007l ; iFeng et al . 2008), and 
used t o probe the turbulent behavior of IGM at hig h red- 
shifts (|Lu et al.|[2O10l ; IZhu. Feng fc Fand[2O10l . l201ll ). The 
last but not the least, searching for redshifted Lya photons 



from star forming galaxies at high redshift is believe d to be 
a bas i c tool to explore the epoch of reionization (e.g. iHaved 
|2010| ; iLehnert et al.|[2uTol ). Therefore, it is crucial to have 
a complete understanding of the radiative transfer of Lya 
photons caused by their resonant scattering with neutral 
hydrogen atoms. 

The radiative transfer of Lya photons in a medium 
consisting of neutral hydrogen atoms has been extensively 
studied either analytically or numerically. Yet, there are 
very few solutio ns on the time-dependent behavior of Lya 
photon transfer (iFieldl Il958l : iRvbicki fc DeriAntoniolll994l ; 
Iffiggins fc MeiksmM 2009). All other analytical solutions are 
time-independent based on the Fokker-Planck e quation 
l|Harringtodll973l ; lNeufeld1ll990l : [Dlikstra et al.ll2006l ). Time- 
independent solutions are important but can only be used to 
describe the "limiting asymptotic behaviors" o f the radia- 
tive transfer (| Adams! Il975l ; iBonilha et ailll979h . They tell 
us nothing about the time scales of the radiative transfer 
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of Lya photons. These time-independent solu tions can not 
descr i be the Wouthuys en-Field (W-F) effect jWouthuvsenl 
Il952l : iFieldl 1 19581 . Il959h . which is essential for the 21 cm 
emission and absorption o f neutral hydrogen at high red- 
shift (e.g. lRov etai]|2009rj ). It is because the Fokker-Planck 
approximation would miss the detailed balance relation of 
resonant scatter ing, which is n ecessary to keep the W-F local 
thermalization !|RvbickJl2006h . 



Numerical method based on the Monte Carlo 
(MC) simulation is also popular i n solving the 
trans fe r of resonant photon s (e. g. lLoeb fc Rvbickl 



1999; Zheng & Miralda-Escudc 2002; Tasitsiomi 


2006; 


Verhamme et alJ 


2006; 


Laurscn & Sommer-Larsen 


2007; 


Diikstra & Loebl 


2008; 


Picrleoni et all 120091: IXu & Wul 


201fj). However, there are verv few works dealing with 



time-dependent problems. For instance, the time-scale of 
the W-F local thermalization is still absent in these studies. 



Time-dependent behavior of the Lya photon transfer 
is especially important to understand observations of Lya 
photons at high redshifts. First, either the life time or the age 
of photon sources at high redshift is generally short. Second, 
the optical depth of the IGM or the halo cloud around the 
sources is generally large. If the time scale of the transfer of 
Lya photons is comparable to the life-time or the age of the 
photon source, the "limiting asymptotic state" will never 
be approached. One more problem is caused by the cosmic 
expansion, of which the time scale is short at high redshifts. 
When the time scale of the cosmic expansion is comparable 
to that of the resonant photon transfer in optically thick 
medium, time- independent solutions probably do not exist. 

Recently, a state-of-the-art numerical method based on 
the WENO scheme has been introduced to solve the integro- 
differenti al equation of th e radiative transfer of re sonant 
photons (|Rov et al.ll2009al lbl; IRov. Shu, fc Fandl201Gi ). It re- 
veals many interesting features of the transfer of Lya pho- 
tons in an optically thick medium, which cannot be seen 
with the time-independent solutions of the Fokker-Planck 
approximation. For instance, it shows that the time scale of 
the formation of the W-F local thermal equilibrium actu- 
ally is short, only about a few hundred times of the reso- 
nant scattering. The double peaked frequency profile of Lya 
photon can not be described by time-independent analyt- 
ical solutions unless the optical depth of v§ photons is as 
large as about 10 6 . This result directly indicates the need of 
time-dependent solutions. 

The WENO algorithm of the integro-differential equa- 
tion of the radiative transfer is fine, but like other high order 
scheme with fixed grid without artificial-viscosity, the com- 
putation time is much more than the Monte Carlo method. 
Therefore, the WENO method would be not easy to deal 
with cases of medium with very high optical depth. The goal 
of this paper is two-fold. The first is to show that the Monte 
Carlo simulation method can properly match the results of 
WENO method on the time-dependent features of moderate 
optical depth. Secondly, we study the time-dependent solu- 
tions of Lya photons escaped from optically thick medium. 
We will not work on specific objects, but focus on the general 
time dependent features which will affect the observability 
of the Lya sources embedded in, or behind optically thick 
medium. 

The paper is organized as follows. Section 2 presents 
the basic models we will study. The Monte Carlo simulation 



method and its tests for time-dependent problems of Lya 
photon transfer in an optically thick medium will be given in 
Section 3. The major results of the time-dependent solutions 
of Lya photon emission in the DLA and IGM models are 
given in §4 and 5, respectively. Problems of absorption by 
optically thick medium are presented in §6. Discussions and 
conclusions are given in §7. All the relevant formulae of the 
radiative transfer of Lya photons and the details of MC 
simulation are presented in Appendix. 



2 PROBLEM 

We study the time dependent transfer of Lya photons in two 
typical models of neutral hydrogen HI medium. The first 
one is the so-called DLA (damped Lya system) halo model, 
in which a source is surrounded by a static spherical halo of 
physical radius r p , consisting of homogeneously distributed 
neutral hydrogen with number density tihi and temperature 
T. The second one is a source at redshift (1 + 2) = 10 lo- 
cated in homogeneously Hubble expanding IGM, of which 
the density and temperature are given by the parameters 
of the concordance ACDM universe. We call it IGM model. 
The radiative transfer equation of the two models are given 
in Appendix A. 

In order to compare with the WENO solutions, for 
the DLA halo model we use dimensionless time and ra- 
dial coordinate, defined, respectively, as r\ = cnmcot and 
r = nHioor-p, where t and r p are the physical variables of 
time and radial coordinate, and oaj^pK is the cross section 
of scattering at the resonant frequency vq = 2.46 x 10 15 s _1 . 
Therefore, r] and r are the time and length in the units of 
mean free flight-time and mean free path of photon vo, re- 
spectively. The value of r actually is equal to the optical 
depth of the spherical halo from r = to r at frequency 
vq. For a signal propagating in the radial direction with the 
speed of light, we have r = r\ + const. 

As usual, in frequency space, we use variable x = 
(y — vo)/Al>d, where Aun = vqVt/c = vo^J 2ksT / m-nc 2 = 
1.06 x 10 n (T/10 4 ) 1/2 Hz is the Doppler broadening at fre- 
quency Vq by thermal motion vt of gas with temperature T. 
The variable x is then the deviation of frequency v from vq 
in units of the Doppler broadening. 

With the dimensionless variables, the specific number 
intensity of photons is I(rj,r,x, /x), where \i = cos# being 
the direction relative to the radial vector r. Thus, all the 
solutions of 1(77, r,x,n) do not refer to a specific density nm 
and size r p (see Appendix § A) . This helps to see the common 
features of the DLA halo model. 

The optical depth of a halo or cloud with column density 
./Vhi at frequency x is 

t(x) = N m v(x) = T <f>(x,a) (1) 

where a(x) is the scattering cross section, tq = Ahicto, and 
4>{x, a) is the normalized Voigt profile given bjQ 

1 Due to different normalization scheme of the Voigt function 
of Eq. (2), our definition of <ro is different from the one used in 
some literatures by a factor 7T 1 / 2 . Consequently, the expressions 
for mean flight time, mean free path and optical depth may be 
different by a factor of 7T 1 / 2 . 
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where the parameter a is the ratio of the natural to 
the Doppler broadening. For the Lya line, a — 4.7 x 
l(n 4 (T/10 4 )~ 1/2 . The profile Eq.(2) describes the joint ef- 
fect of the Gaussian distribution of the velocity of neutral 
hydrogen atom and the Lorentz profile of cross section in 
the rest frame of the atom. For an expanding (or collaps- 
ing) halo or turbulent gas cloud, the bulk velocity of the gas 
might be larger than the thermal velocity vt- Even in this 
case, the Doppler or thermal broadening is still important, 
as it is the key factor leading to a local thermal equilibrium 
of Lya photons (the W-F effect). 

In an optically thick spherical cloud, most time- 
dependent behaviors can be described by two time- 
dependent distribution functions: 1.) angularly aver- 
aged, r 2 -rescaled, specific number intensity j(rj,r,x) = 
(r 2 /2) f. I(r),r,x, tl)dfx, and 2.) r 2 -rescaled number flux 

f(r/,r,x) = (r 2 /2) [il(r),r,x, fj,)dfj,, which describes the 
photons escaped from a spherical halo of radius r. The equa- 
tions of j and / are given in Appendix Eqs.(Al) and (A2). 

For the IGM model, we use the physical variable t and 
r p because of the specific cosmological model we adopted. 
The temperature of IGM is taken as T = 100 K and the 
parameter 7 = 1/tqp [see Eq.(A3)], which describes the 
Hubble expanding, is equal to 1.4 x 10 -6 in the concordance 
ACDM universe. 



3 METHOD AND TEST 

3.1 Time dependent Monte Carlo simulations 

We use Monte Carlo (MC) method to simulate j(r/, r, x) and 
f(rj, r, x) of the previous section. Most Monte Carlo codes of 
simulating time-independent solutions of radiative equation 
can easily be modified to deal with time-dependent prob- 
lems. If the feedback of photon transfer on the parameters 
of neutral hydrogen is negligible, the radiative transfer equa- 
tions are linear with respect to I, and thus to j and /. 
Thus, linear superposition of the solutions of the sources 
is valid. The time-dependent solutions can simply be given 
by a weighted summation over the results of a single flash. 
The weight of the summation is proportional to the time- 
dependent flux of the light source. 

We w i ll em ploy the same MC algorithm as used in 
IXu fc Wul |20ld ). of which some of the details are given in 
Appendix B and C. The major modification from th e ear- 
lier methods, such as lZheng fc Miralda-Escudd |2002), is to 
record the time at each collision of the photon, and to take 
snapshots of photon distribution in spatial and frequency 
space based on these time stamps. Thus, the time-dependent 
solutions j and / of an arbitrary source can be given by a 
synthesis of the fluxes at different epochs from a single flash 
source. 



3.2 Test with ^-profile of flux 

As a test of the MC method, we calculate the flux f(rj, r, x) 
at the outer boundary ro = 10 2 of a DLA halo. This result 
is shown in Fig. 1, in which the curves are the MC results 
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Figure 1. The Lya photon flux f(r],r,x) from a DLA halo with 
radius ro = 10 2 and parameter a = 10 — 3 . The curves are from 
MC simulations at time rj = 500, 2000 and 3000, respectively 
from bottom up. The last two curves (2000 and 3000) are already 
overlapped with each other. The symbol points marked by "+" 
and "x" are from t he WENO numerical so lution of the radiative 
transfer equation l|Rov. Shu, fc Fandlioioh at epochs 77 = 500, 
2000, respectively. 

with parameter r\ = 500, 2000 and 3000. They show typical 
double-peaked profile. The curves of 77 = 2000 and 3000 ac- 
tually are the same. The overlapping indicates that f{r], r, x) 
is already at the limiting asymptotic state, or saturated for 
time 77 ^ 2000. In Fig. 1, the data points are given b y 
the WENO numerical solutions of lRov. Shu, fc Fand l|2010h . 
Therefore, the MC method can match the time-dependent 
WENO solutions. 

In the WENO method, the solutions of the angularly 
averaged specific intensity j(rj,r,x) and flux f(rj,r,x) at 
boundary ru should satisfy the condition j(-q,ru,x) — 
2/(/7, ro, x) (Unno 1955). This is the result of Eddington 
approximation and the assumption of no incoming photons 
at the boundary. Our MC simulation results of jr (77, ro, x) at 
time 77 = 500 are shown at various radius in Fig. 2. At the 
surface where r = 100, we see j « 2/ at the center frequen- 
cies when compared with Fig.l. It shows that the MC sim- 
ulation can pass the test of Unno's boundary condition. We 
also find that j — 2/ relation is only valid at the boundary 
at center frequencies. Slightly beneath the surface, at radius 
r = 99, 98, we find that j ~ 4/ , and 5/ , respectively. The 
enhanced photon intensity is a result of backward scattered 
photons near the boundary. 

Besides the curve of r = ro — 100, Fig. 2 also plots 
j(v,r,x) at time 77 = 500, but for r = 80, 90, 95, 98, 
99. We see that all the curves of r < 100 are almost flat 
in the range of \x\ < 2. It means that the frequency dis- 
tribution of photons is thermalized near the resonant fre- 
quency vq- That is, the frequency distribution is of Boltz- 
mann tx exp(— h/ ^P x), where T is the kinetic temperature 
of neutral hydr ogen gas in the halo. This is the W-F local 
thermalization (|Wouthuvsenlll952l ; iFieldl Il958l . 119591 1. Fig. 
2 tells us that the W-F local thermalization is achieved by 
resonant scattering even at r = 99. Yet photons at the out- 
ermost layer r — ro have not yet been thermalized, as the 
optical thin layer does not carry enough number of scatter- 
ing. This result is similar to the WENO solution (|Rov et al.l 
l2009bl : iRov. Shu, fc Fan3l2010h . 
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Figure 2. The MC simulation results of the angularly averaged 
specific intensity j(r), r, x) at time rj = 500 for DLA halo model of 
parameter a = 0.001 with size rpj = 10 2 . The curves correspond 
to r = 100, 99, 98, 95, 90, 80, respectively from bottom up. The 
flux / at time rj = 500 at the boundary r = 100 from Fig.l is 
marked with " +" symbols for comparison. 



3.3 Test with [x±[-ro relation 

The second test is on the a;±|-ro relation, where x± are fre- 
quencies of the two peaks of the double-peaked profile as 
shown in Fig. 1, and ro is the optical depth parameter de- 
fined in Eq.(l). This relation has been studied by many time- 
independent solutions based on the Fokker-Planck equa- 
tion. The major conclusion is the so-called (aro) 1 ^ 3 - law, i.e. 
x± = ±A(aTo) 1 / 3 , where A is a constant of order 1 | Adams 
1971 1 19751 ; iHarringtonl Il973l ; iNeufeldl fl99fj| ; iDiikstra et al 
200£ ). It is well known that the (ar) 1 ^ 3 -law is available only 
when the optical depth ro is very large. However, the x±- 
to relation available for various ro has not been calculated 
until very recently. It may be due to the absence of proper 
numerical solver of the integro-differential equation of reso- 
nant photon transfer in optical thick medium. The WENO 
solver provided the first x±-to relation of DLA halo with 
moderate and high optical depth. 

The x±-tq relation given by the MC simulation is pre- 
sented in Fig. 3, which shows that the (ar) 1 ^ 3 -law is sig- 
nificantly deviating from the MC results when aro ^ 10 2 , 



Figure 3. The two-peak positions x± vs. (aro) given by MC 
simulation (data points) of a DLA halo with optical depth ro. 
The parameter a is taken to be 5 X 10 — 4 (bold line) and 5 X 10~ 3 
(light line). The dashed straight line of log x±-log aro with slope 
1/3 is to show the (aro) 1 / 3 -law. 



where the optical depth is smaller than to — 2 x 10 at 
a = 5 x 10~ 3 . This result is the same as the WENO solu- 
tion. In Fig. 3 the parameter ra nge of [aro] is larger tha n 
that of WENO solution [Fig. 4 of lRov, Shu, fe Fane) (|2010T l]. 

In the range 10 -2 < aro < 10 2 , the |a;±|-ro relation is al- 
most flat with \x±\ ~ 2. It is because the double-peaked pro- 
file is from photons stored in the frequency range of \x\ < 2 
and in the local thermal equilibrium state. The positions 
of the two peaks, x±, actually is about the same as the 
frequency range of the local thermalization. The frequency 
width | a; | ^ 2 of the local thermal equilibrium state is deter- 
mined by the Doppler broadening, and very weakly depen- 
dent on ro. Thus, once the photons in local thermal equi- 
librium state are dominant, we always have x± ~ ±2. This 
point can also be seen with Figs. 1 and 2, in which the po- 
sitions of the two peaks are kept to be \x±\ ~ 2 despite that 
the intensity of the flux increases with time significantly. 

When aro < 10" 2 , the curves of |x±| are no longer 



determined by one variable aro, instead by variables a and 
ro, separately. The value of \x±\ shows a quick drop to zero 
at aro ~ 10~ 2 for a = 5 x 10~ 3 , and aro ~ 10~ 3 for a = 5 x 
10~ 4 . The W-F thermal equilibrium cannot be established 
in halos with small optical depth ro << 10 2 , and therefore, 
photons from these halos do not have double-peaked profile. 

Since thermalization will erase all frequency features 
within the range \x\ ^ 2, the double-peaked structure 
doesn't retain information of the photon frequency distribu- 
tion within \x\ < 2 at the source. It is impossible to probe 
the frequency profile for \x\ < 2 Lya photons of the source 
from the escaped Lya photons. This property can also be 
used as a test of simulation code. That is, the simulation 
results should be independent of the profile of Lya emission 
from the sources, only if the profile is non-zero within the 
range \x\ < 2, i.e. it should not matter whether the source 
is monochromatic, or has a finite width around vo- 



4 LYa PHOTON EMISSION: DLA MODEL 
4.1 Time dependent Lya escape 

It is well known that the spatial transfer of Lya photon in 
optically thick halo is not simply a Brownian random walk. 
The time scale of Lya photon escape from optically thick 
halo is much shorter than that of Brownian diffusion. It is 
because the spatial transfer depends on the diffusion in fre- 
quency s pace. This is the so-called "single longest excursion" 
process if Adamsll 19721) . However, earlier estimates of the es- 
caping time scale based on "single longest excursion" can 
not describe the details of the time-dependent behavior of 
photons escaping from optically thick medium. 

If the central source of a DLA halo is assumed to be a 
photon flash, the time dependence of the luminosity of the 
photon source is proportional to S(rj). Without scattering, 
the luminosity at the boundary of the halo with size r should 
still be a delta function as oc <5(jy — r), i.e. it is also a flash, 
but with a retarded time r, which is the time needed for a 
freely streaming photon from the center to the edge of the 
halo with speed of light. Considering the effect of resonant 
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Figure 4. The light curve of a flash of Lya photons in a DLA 
halo model of a = 0.0005 with optical depth to = Hd = 2 X 10 7 . 
One photon is supposed to be emitted from a central source at 
time r\ = 0, i.e. the integral of the light curve is equal to 1. 



scattering, the luminosity of escaped photons will no longer 
be a flash. The light curve of the luminosity from such a 
source for a halo with size r = 2 x 10 7 and a = 5 x 10~ 4 
is shown in Fig. 4, in which F(rf) = j f(r],x)dx is the flux 
integrated over all frequencies of escaped photons. 

Fig. 4 shows that the light curve lasts from time r\ ~ f 8 
to 5 x 10 s . The peak of light curve F max is at 77 max ~ 2 x 10 8 , 
and the time duration Ar/ of F{rf) > F max /2, is also about 
2 x fO 8 . Since the source of photons does not contain any 
time scale, both the numbers ?7 max and A77 are from the size 
or optical depth of the halo. The time scale ?7 max means that 
the retarded time of the escaping of photons from the halo 
is as large as ~ lOr or 10ro. The amount Ar/ means that the 
time-distribution of the escaped photons are significantly 
spread out from a Delta function 5(77) to Ar] ~ fOr. 

Without resonant scattering, photons emitted from the 
source will escape from the halo at time r\ = r = to. With 
resonant scattering, the majority of photons emitted from 
the source will not escape from the halo until time r\ = 
?7max ~ 10 s . Therefore, a huge number of resonant photons 
is stored in the halo. The resonant nature of Lya photon 
scattering let the photons to stay in the halo with a time 
scale equal to 10 times of the optical depth to of the halo. 

In our model the destruction processes of Lya photons, 
such as the two-photon process, are not considered, and dust 
absorption is ignored too. The number of photons is con- 
served. Thus, we have J F(rf)dr) = 1. Therefore, the light 
curve, F{rf) of Fig. 4 can be understood as the probabil- 
ity distribution of the time of Lya photon escaped from a 
r — 2 x 10 7 halo. With dimensionless variables, the curve of 
Fig. 4 is also the probability distribution of the total length 
of the path of a photon transferring from r = to r = 2xl0 7 . 
In this context, 7? max and Ar/ can be used as the most prob- 
able length of the path, and the variance of the distribution, 
respectively. Considering that the most probable path length 
Wmax and the variance Ar) have the same order, the spatial 
transfer of Lya photon in optically thick halo essentially is 
still a random process of diffusion. 

The light curve of a flash source (Fig. 4) can easily be 
generalized to arbitrary sources. Considering the total flux 
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Figure 5. Time dependent solutions of escaping coefficient 
/osc(f?) for photon sources located in the center of a DLA halo 
of T = 10 4 K with optical depth r = 2 X 10 7 . The life dura- 
tion of the source ?7iif c is taken to be ?7ijf e = 10 6 , 10 7 , 10 8 , 10 9 
and 10 10 , respectively, from bottom up using color code orange, 
violet, cyan, green, and black. 

is linearly dependent on sources, the total flux L(rf) is given 
by 

L(n)= f V ' Fin-n'^dr,' (3) 

Jo 

where F(rj) is the curve of Fig. 4, and s(rj) is the time- 
dependent Lya photon flux of the sources. 

4.2 Escape coefficient 

For a source with stable luminosity of Lya photons, escape 
coefficient of Lya photons emergent from a halo with size r 
can be calculated by f eBC (r],r) = F(rj)/F = J f(rj,r,x)dx if 
the luminosity of sources is normalized. Since the number of 
Lya photons is conserved, the total number of escaped pho- 
tons should be equal to the total number of emitted photons 
by the source when a stable state is reached. Thus, the es- 
cape coefficient of a stable source should reach to 1 when 
r\ is large enough. However, before the system approaches 
to stable state, the escape coefficient / eS c (??, r) can be much 
less than one. 

We have calculated the time dependent solution of the 
escape coefficient f e sc(rj) f° r a stable photon source with 
unlimited life time located at the center of the halo with 
optical depth to = 2 x 10 7 . At time r; max = 2 x 10 s , the 
escape coefficient f esc (ri) ~ 0.37. ff we define the time scale 
?7 3a t of reaching saturated or stable state as / eS c(^ sa t) = 0.95, 
it yields r/ Ba ,t = 6.1 x 10 8 , or ~ 30ro. Therefore, / csc is always 
significantly less than 1, when the age is less than about 30to 
or 30 r. 

Fig. 5 presents the time dependent solutions of escape 
coefficient / e8C (??) for photon sources in a T = 10 4 K neutral 
medium with limited life time to be 7/ufe = 10 6 , 10 7 , 10 8 , 10 9 
and 10 10 , respectively, from bottom up. ft shows that the 
escape coefficient is always less than 1 if rjute < 10 7 . In 
this case, the time scale of the escaping of Lya photons 
in ro = 2 x 10 7 halo is much larger than the life time of 
the source, and therefore, photons emitted within the time 
duration r\ = to rj = r/ut,, are fully spread over a time scale 
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Figure 6. Time dependent surface brightness of a central source 
surrounded by a DLA halo of T = 10 4 K with optical depth to = 
2 X 10 7 . The source emits Lyct photons with stable flux L = 10 42 
erg s — 1 in the time range from t = to t = 10 6 years. The size 
of DLA halo is taken to be rrjLA = lkpc, 10 kpc, 100 kpc and 
1 Mpc, respectively from top down using color codes black, red, 
orange and green. 



Arj 3> r?iif c - Thus, the escape coefficient of a source with 
short life time is always much less than 1 when the DLA 
halo is optically thick. This mechanism would be important 
for understanding the observable features of high redshift 
objects such as GRBs or first stars. Even when dust absorp- 
tion is negligible, the escape coefficient can be small, even 
very small, when the age of the object is small. 

It is interesting to see that although different curves of 
Fig. 5 correspond to very different rjnf e , all the curves with 
7/iif c < 10 7 approach their maximum at about the same time 
r\ ~ 2x 10 s . This time scale actually is the 77 ma x shown in Fig. 
4. The stored photons yield a delayed emission with time 
scale of about 77 max ~ 2 x 10 8 . The delay time is independent 
of the life-time Tjnfe of the source, but dependent only on the 
optical depth of the halo. 



4.3 Surface brightness and the size of DLA halo 

For a stable source, the time-independent surface bright- 
ness (SB) of a stable source located in DLA halo is simply 
proportional to r^ A , where tdla is the radius of the halo. 
When the life time of the source is short, the time-dependent 
behavior of the surface brightness, SB(t), is not simply pro- 
portional to r~^LA- To demonstrate this point, we calculate 
the surface brightness of a source with life time 10 6 years and 
Lycv luminosity L — 10 42 erg s _1 in a medium of a = 0.0005. 
The results are presented in Fig. 6, in which the size of the 
halos is taken to be tdla = lkpc, 10 kpc, 100 kpc, and 
1 Mpc. In Fig. 6, the optical depth of the DLA halo sur- 
rounding the source is assumed to be to = 2 x 10 7 , which 
correspondes to a typical DLA halo with column density 
JVhi = 2 x 10 20 cm" 2 (Eq.(l)). 

From Fig. 6, we first see that the curves do not show the 
behavior of being proportional to r^ A . The maximum sur- 
face brightness of halo with tdla = 1 kpc is about 10 -17 erg 
s _1 cm -2 arcsec -2 , while the maximum surface brightness 
of tdla = 1 Mpc is only ~ 10 -25 erg s -1 cm -2 arcsec -2 . 



That is, when r~, 2 A decreases by a factor 10 6 , the maximum 
surface brightness descreases by a factor 10 s . 

Secondly, the shapes of the curves of Fig. 6 for different 
»"dla are very different. The curve of tdla =1 kpc shows 
saturation at t > 3 x 10 4 years, while all others don't have 
a saturated phase. 

These time-dependent behaviors are also due to the 
time scale of Lycv photon transfer. For a halo with given 
column density Nm, or to, the time scale of the Lya pho- 
ton transfer is about r) max ~ 10ro, which is independent 
of the size tola of the halo. However, r; m ax ~ 10ro corre- 
sponds to a physics time i max = TTmaxTDLA/c, which is r D LA- 
dependent. For DLA halo with tdla =1 kpc, the time scale 
tmax = »7maxfDLA/c ~ 3 x 10 4 years. It is much less than 
10 6 years, and therefore, the surface brightness is saturated 
when t>3x 10 4 year. For DLA halo with tdla =100 kpc, 
the time scale t max = ^max^DLA/c ~ 3x 10 8 years. It is larger 
than the life-time of the source, and therefore, the surface 
brightness can not approach a saturated state. For DLA halo 
with r D LA > 100 kpc, the time scale t max = JfaaxT'DLA/c 
would be much larger than 10 6 years. The source is more 
like a single flash and its radiative transfer is highly time 
dependent (Fig. 4). 

Thus, we may conclude that for DLA halos with a given 
column density, the surface brightness will be proportional 
to r^LA an d K > 2. This result would be useful to estimate 
the size of the DLA halo with observed surface brightness 
and the model of the sources. 



5 LYa PHOTON EMISSION: IGM MODEL 

5.1 Lycv escape 

The mechanism of the escaping of Lya photons from ex- 
panding opaque IGM at high redshift is different from that of 
DLA halos. For the former, besides the diffusion in frequency 
space caused by the the resonant scattering, we should also 
consider the frequency redshift caused by cosmic expansion. 
Photons will escape from Gunn-Peterson trough, once their 
frequency is redshifted enough. 

Similar to the previous section, we calculate the es- 
cape coefficient using the number of escaped photons. Fig. 
7 presents the time-dependent solutions of the escape coeffi- 
cient /esc(i) for sources at redshift 1 + 2 = 10 surrounded by 
the IGM of Hubble expanding universe. The number density 
of neutral hydrogen is given by the cosmology parameters 
of the ACDM model. The IGM is assumed to be not reion- 
ized yet. The time duration tu{ e of photon source has been 
explored with values t ufc = 10 3 , 10 4 , 10 5 , 10 6 , 10 7 , and 10 8 
years. For sources at 1 + z = 10, tm e cannot be larger than 
the age of the universe ~ 10 8 years. The source is starting 
to emit photons at t — 0. 

The curves of Fig. 7 look very similar to those of Fig. 5. 
Therefore, one can explain Fig. 7 in the same way as Fig. 5, 
replacing the optical depth of the DLA halo with the Gunn- 
Peterson optical depth (Eq.(A5)) of the expanding IGM at 
(I + 2) = 10. In Fig. 7, all the short life time curves with 
tiifc < 10 7 years reach their maximum at about the same 
time t max ~ 10 6 years, which is the time scale of Lya pho- 
ton escape from the Gunn-Peterson trough in an expanding 
IGM at 1 + z = 10. One can then conclude that the escape 
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Figure 7. Time dependent solutions of escaping coefficient f eBC (t) 
for a source at redshift 1 + z = 10 surrounded by the IGM of 
T = 100K of the expanding universe. The time duration of the 
emission of the source is taken to be tvif e = 10 3 , 10 4 , 10 5 , 10 6 , 
10 7 , and 10 8 years, respectively from bottom up using color codes 
blue, cyan, green, yellow, orange, and black, t is the delayed time 
of an escaped Lya photon with respect to the first escaped free 
streaming photon from the source, measured at the redshift of 
the source. 

coefficient of sources with life time less than t max ~ 10 6 years 
should always be less than 1, even without dust absorption. 

The curve of the shortest life time (inf e = 10 3 ) of Fig. 7 
can be thought to represent the light curve of a flash source 
in IGM model, like that of Fig. 4 for DLA model. The curve 
of tnfe = 10 3 of Fig. 7 has a long tail. The long tail basically 
is a power law of t, and is a joint result of radiative transfer 
and Hubble expansion velocity field. As mentioned in §4.1, 
for the light curve of Fig. 4, the maximum ry ma x and variance 
(or the width of the light curve) A77 are about the same. Yet, 
the power law long tail leads to the maximum t max of the 
curve of tiif e = 10 3 of Fig. 7 to be much less than the width 
of the curve. Since the suppression of escape coefficient is 
mainly given by the width of the light curve. Therefore, the 
power law long tail of Fig. 7 implies that the stable state of 
/esc = 1 takes a much longer time to approach, or can never 
be approached within the age of the universe. 

Although Fig. 7 shows that the escape coefficient of 
sources with tijf e ^ 10 6 years is larger than those of sources 
with tijfe ^ 10 6 , it does not mean that the former is easier to 
be observed than the later. This point can more clearly be 
revealed with the time-dependent solution of surface bright- 
ness. 



5.2 Surface brightness 

We calculate the mean surface brightness SB(t) = \ — — 

defined as the flux divided by the averaged area of the 
source, where r esc is the projected distance to the source 
on the sky at which a photon escapes. Thus, ty < r 2 sc >(t) 
is the mean of the area at time t on the plane perpendicular 
to the line of sight. The result is presented in Fig. 8, which 
plots the SB(t) for sources with Lya luminosity (erg/sec) 
and life time (year) paired as: (10 4e , 10 3 ), (10 45 , 10 4 ), (10 44 , 
10 5 ), (10 43 , 10 6 ), (10 42 , 10 7 ), and (10 41 , 10 8 ). That is, all 
sources emit the same amount of 2 x 10 67 Lya photons in 
their whole life span. The sources start to emit at time t = 0. 
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Figure 8. Time dependent surface brightness from a central 
source surrounded by IGM of T = 100 K in expanding universe. 
The source is assumed to be located at redshift 1 + z = 10, and 
emitting a total of 2 X 10 67 Lya photons from time t = 0. The 
time duration of the emission of the source is taken to be ii;f c = 
10 4 , 10 5 , 10 6 , 10 7 , and 10 8 years, respectively from top down using 
color codes cyan, green, yellow, orange, and red. The flash source 
(*lifc = 10 3 years) is drawn with black dashed line. 



Fig. 8 has two remarkable features. First, the four curves 
of (10 46 , 10 3 ), (10 45 , 10 4 ), (10 44 , 10 5 ), and (10 43 , 10 6 ) are 
almost the same. The common curve has a maximum at 
t ~ 10 6 year and then starts decaying with a power law 
SB(t) oc t~ 4 . The time scales of the maximum and the 
power law tail are about the same as that of the light curve 
of Fig. 7 with ime = 10 3 . Therefore, for sources with life 
time less than t ~ 10 6 years, the surface brightness is only 
dependent on the total number N of Lya photons emitted 
from the source, regardless of their life time. This is be- 
cause the photons should wait for about t ~ 10 6 years before 
their escaping from the IGM. The stored photons are locally 
thermalized, and therefore, the information of the "age" of 
the photons emitted at time < 10 6 will be forgotten during 
the thermolization. This property actually is also valid for 
the sources of (10 42 , 10 7 ), and (10 41 , 10 8 ). For these two 
cases, the total numbers of Lya photon emitted in t ^ 10 6 
years are smaller, respectively are 2 x 10 66 , and 2 x 10 65 
and therefore, their SB(t) at t ~ 10 6 is less than that of 
N — 10 67 by factors 10 and 100, respectively in the figure. 
The curves of tuf e = 10 4 , 10 5 years become scalable from 
~ 1 Myr. The tuf e = 10 6 yr curve is just starting to move 
away from the flash source solutions. The relaxation time 
looks to be around 0.3 — lMyr when the flash source reaches 
its maximum. This is in agree ment with the estimate from 
iRvbicki fc Dell'Antoniol 1 19941 ) where the relaxation time is 
(-%)* * ts ~ 0.3 Myr for our chosen redshift (1 + z = 10) 
and temperature (T = 100K). Fig. 8 shows various scaling 
relations. For example, for the flash source of tuf e = 10 3 yr, 
the asymptotic slope sets in from ~ 1 Myr as a result of ra- 
diative transfer in Hubble expansion velocity field. For the 
tufe = 10 8 yr curve in the figure, another straight part sets 
in from 1 Myr to 50 Myr even before the asymptotic slope is 
reached, as a joint result of photon emission, photon transfer 
and Hubble expansion. These time scales are comparable to 
the possible lif e time of the sourc e s, wh ich was previously 
pointed out bv lHiggins fc Meiksinl (|2009l V 
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The second feature of Fig. 8 is the monotonous decrease 
of SB(t) when t is large, regardless of the life time of the 
source. It is very different from all the curves of Fig. 5, which 
will approach a stable or saturated state if the life time of 
the source is long enough. For any sources in expanding 
neutral IGM, the mean surface brightness SB{t) does not 
approach a saturated or stable state. In other words, a time- 
independent solution of the surface brightness doesn't exist 
in the IGM model. This is simply because the increase of 
volume of the expanding IGM, which stores Lya photons, 
is faster than the number of Lya photons redshifted to fre- 
quency x < — 2. This feature is similar to the evolution 
of ionized halos around a UV photon source in expanding 
universe. The ionized radius can never approach a stable 
state required by a Stromgren sphere, because the increase 
of the i onized radius is always lower than the comoving 
velocitv (|Shapiro fe Girouxlll987r ). Therefore, in Fig. 8 there 
is no flat section for every curve. The maximum of the sur- 
face brightness of a Lya photon source with luminosity L at 
1 + z — 10 scales approximately as ~ 10 -21 x (L/10 43 erg/s) 
erg s _1 cm -2 arcsec. We should emphasize that the max- 
imum can be reached only for sources with age equal to 
about 10 6 years. The surface brightness would be less than 
the maximum when the source age is younger or older than 
10 6 years. That is, in term of surface brightness, a stable 
source will yield a time-dependent curve. 
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Figure 9. Time dependent red damping wings formed in a DLA 
halo of parameter a = 0.0005 with size r and optical depth to set 
to be r = ro = 2 X 10 7 . The central source starts to emit photons 
with a flat spectrum at r) = 0. The profiles of the red damping 
wing are shown at times r} = (1 + y)Vai where r\ c = r = ro is the 
time of free flight from r = to r, and y is taken to be 0.01, 0.5, 
2, 5, 10, and 50, respectively from bottom up using color code 
violet, blue, orange, yellow, green, and cyan. The black curve, 
which is almost the same as the violet curve, is given by the 
Voigt absorption, i.e. f(x) = e~ T ^ and t(x) is from Eq.(l). 



6 LYq photon absorption 

6.1 Time dependent red damping wing: DLA 
model 

A typical problem of Lya absorption at high redshift is the 
red damping wing, or HI damping wing of high redshift 
sources. The optical aft erglow of GRBs at high redshifts has 
show n this feature(e.g. iTotani et al.|[200r3 ; ISalvaterra et al.l 
2009). We calculate the time-dependence of the red damp- 
ing wing of a DLA halo with radius and optical depth ro to 
be r — to = 2 x 10 7 . The frequency spectrum of the central 
photon source is assumed to be flat with flux equal to one. 
If the source starts to emit photons at time r\ = 0, photons 
without undergoing scattering or collisions will escape from 
the halo at time r\ c = 2 x 10 7 . The profiles of red damping 
wing at times later than r\ c are plotted in Fig. 9, in which 
we take rj = (1 + y)Vc, where y = 0.01, 0.5, 2, 5, 10 and 50. 

Fig. 9 shows that the red damping wing at time 77 ^ 
l.Olr/c can be well described by a Voigt profile f(x) = e~ T< - x \ 
where t(x) is given by Eq.(l). Therefore, it would be fully 
reasonable to fit the red damping wing of GRB's optical 
afterglow with a Voigt absorption profile, because the red 
damping wing is measured only a few hours or a few days 
after the GRB explosion, the time 77 is very close to r] c . Even 
when 77 = 1.5rj c (the blue line in Fig. 9), the red damping 
wing can still be approximately fitted by a Voigt profile. 
However, the fitting at time 77 — 1.5rj c will yield a smaller 
ro than that of the fitting at 77 = 1.01?7c. Therefore, the 
column density of HI atoms given by the fitting at 77 = 1.5?7 C 
is underestimated. 

The red damping wing at 77 5r/ c shows a shoulder with 
frequency similar to that given by Fig. 3 with parameters 
a = 5 x 10~ 4 and ro = 2 x 10 7 . Therefore, the shoulder 



is from the stored resonant photons. Fig. 9 shows that the 
curve of red damping wing has become saturated at time 
77 ^ 5077 c . This is larger than the time scale ~ 10?7 C shown 
in Fig. 4. It is because the photons with frequency |x| > 2 
of the continuous spectrum can also be stored. According to 
the redistribution function ( A6) , the probability of transfer- 
ring a \x\ > 2 photon to \x\ < 2 is larger than that from 
I a; I < 2 to \x\ > 2, and therefore, in frequency space the 
net effect of the resonant scattering is to bring photons of 
the continuous spectrum background to become Lya. More 
photons be stored leads to the larger time scale of the satu- 
ration. 



6.2 Time dependence of red damping wing: IGM 
model 

The problem of the red damping wing for the IGM model is 
very different from that of the DLA model. If the frequency 
spectrum of source is flat, cosmic redshift will continuously 
provide Lya photons from blue side of the spectrum. Since 
the time scale of the W-F thermalization is much shorter 
than the time scale of the cosmic expansion, the photons 
moved into the frequency range ~ uq from > vq is quickly 
thermalized. Thus, a huge number of th ermalized photons 
is stored in the frequen cy range \x\ < 2 (|Rov et al.ll2009bl ; 
iRov. Shu, fc R ml 1201011 ." This process marks the difference 
between the IGM and the DLA models, and also the dif- 
ference between our IGM mo del and IGM models of some 
others (lLoeb fc Rvbickil I1999T ), in which the source emits 
only uq photons. 

Fig. 10 shows the red damping wings of a source in 
T=100 K IGM at 1 + z = 10 with age t = 0, 10 5 , 10 6 , 
5 x 10 B , and 10 7 years. When the time is less than 10 5 year, 
the Lya photons have not been significantly stored yet, and 
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Figure 10. The time dependence of the red damping wing given 
by simulations in an expanding IGM model. The central source 
at 1 + z = 10 starts to emit photons with flat spectrum at t = 
into IGM of T=100K. The curves correspond to t = 10 5 , 10 6 , 
5 X 10 6 , and 10 7 years, respectively from bottom up using color 
codes blue, orange, cyan, and black. The red line is given by pure 
absorption f(x) = e~ T ^ where t(x) is from Eq.(Cl). 

the profile of red damping wing is about th e same as a pure 
absorp tion, which has been addressed by iMiralda-Escudd 
(1998). Once the time is larger than 10 6 year, the stored 
photons yield a shoulder at peak frequency —x ~ 500. Fig. 
10 also shows that the profile of the red damping wing seems 
not yet to approach a saturated or time-independent state 
at t = 10 7 years. 



turbulence of the IGM field, and the observational aper- 
ture, etc. However, one point seems to be clear that under 
these circumstances, stable solutions will not exist. All these 
effects should be studied with time-dependent solutions of 
the Lya photon transfer. Since some high redshift Lya emit- 
ters have already been observed, static solutions may not be 
enough to understand these observations. 

For the DLA model, the saturated state or the time- 
dependent state, does exist. Therefore, it is reasonable to 
find this state directly by time-independent solution. How- 
ever, time-independent solutions can not tell us the suppres- 
sion of the escape coefficient in the time range of r\ < IOto, 
to being the optical depth of the halo. When to is large, the 
escape coefficient will be substantially suppressed for time 
range comparable to the age of the universe. This is equal 
to say that the escape coefficient is much less than one in 
all time. 

Low escape coefficient is often explained by dust absorp- 
tion. Dust extinction is effective only when the medium is 
optically thick. Therefore, once we use the model of dust ab- 
sorption at high redshift, the time-dependent behavior must 
be considered. A WENO time-dependent solution of the Lya 
photon transfer in dusty medium will be reported soon. 

We didn't consider the Lya photons from the recombi- 
nation in the ionized halo around the center source. Simply 
changing the luminosity of the central light source by adding 
the photons from stable recombination in the halo is a poor 
approximation, as the time scale of the recombination is 
not small. Therefore, with the DLA model, time-dependent 
solutions are also essential for understanding Lya-related 
phenomenon. 



7 DISCUSSIONS AND CONCLUSIONS 

The time-dependent behaviors in physics and frequency 
spaces of the Lya photons emergent from optically thick 
medium have been extensively studied with Monte Carlo 
method. The first conclusion is that time-dependent solu- 
tions are essential not only for strong time dependent sources 
like GRBs, but also for stable sources, especially when the 
cosmic expansion needs to be considered. Cosmic expansion 
makes the radiation transfer equation not invariant with re- 
spect to the transformation of time shift. Consequently, the 
radiative transfer equation doesn't have time-independent 
solutions in principle. Time-dependent behavior are essen- 
tial. 

With the model of IGM in the Hubble expanding uni- 
verse, we show that the time-dependence of the brightness 
of a stable source at high redshift is like that of the light 
curve of an "explosion", i.e. in the first phase, the surface 
brightness is increasing and approaching a maximum, and 
then, in the second phase, it decays monotonously. Although 
our calculation is only based on a model source at redshift 
(1 + z) = 10, we believe that the feature of the nonexistence 
of a time-independent solution, or the nonexistence of a sat- 
urated state of the surface brightness by IGM scattering, 
would be hold for various high redshift sources. 

The IGM model may be too simple and ideal. Many 
effects have not been considered, such as the effects given 
by the inhomogeneity of the HI density distribution, the 
nonuniform ionization, the bulk velocity, the wind and the 
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APPENDIX A: RADIATIVE TRANSFER OF 
LYa PHOTONS 

With dimensionless variables, the specific intensity I of pho- 
ton number is a function of r\, r, x and \i. When the optical 
depth is large, we can take the Eddington appr oximation. 
Radiative transfer equations are (|Rov et al.ll2009cfi 

% + % = ~4>{x\ a)j + f TZ(x, a)jdx' + 7 J- + r 2 S, 
Of] or J ox 

(Al) 

df ldj 2 j df ,, ,.„, 

where j, f are the rescaled quantities of physical inten- 
sity, j{v, r , x ) = r 2 j P = t 2 \ f*i I{v,r,x,fj,)dfJ, is tne an- 
gularly averaged specific intensity and f(rj,r,x) = r 2 f p — 
r 2 | /j,I(ri,r,x, [i)d[j, is flux. The mean intensity j(r),r,x) 
describes the a;— photons trapped in the halo by the resonant 
scattering, while the flux f(jj, r, x) describes the photons in 
transit. The parameter 7 = 1/tgp can be calculated as 
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The term S describes a constant physical photon source. 
Since Eqs.(Al) and (A2) are linear with respect to j and 
/, the solutions of j and / for sources S(r/, x) can be 
given by a superposition of solutions j(rj, r, x; r/o, xo) and 
f(rj,r,x;rio,xo), which are the solutions of the source 5* = 

S(v - rio) s ( x - x o)- 

The resonant scattering is described by the redistribu- 
tion function TZ(x,x';a) which is the probability of a pho- 
ton absorbed at the frequency x 1 , and re-emitted at the fre- 
quency x. It depends on th e details of the scattering (|Henvevl 
Il94ll ; lHummerlfl962l . 1 19691 ). If we consider coherent scatter- 
ing without recoil, the redistribution function with the Voigt 
profile is 



1Z(x, x'; a) 

[' 



tan 



1 

^372 

11 



-x'|/2 



(A4) 



+ u 



— tan 



'-)} 



du 



where x m - ln — min(a;,a; ) and £ max = m&x(x,x ). This re- 
distribution function is normalized as lZ(x, x')dx = 

<f)(x,0) =ir- 1/2 e- x2 . 

Thus, the frequency of photons will be changed during 
the evolution governed by Eqs.(Al) and (A2), while the to- 
tal number of photons is conserved. That is, the destruction 
processes of Lyq photons, s uch as the two-p hoton process 
l|Spitzer fc Greensteinlll95ll ; IOsterbrocklll962l ) and dust ab- 
sorption, are ignored in equations (A1),(A2). 



APPENDIX B: MONTE CARLO SOLVER 

Monte Carlo simulation starts with releasing a photon at the 
source placed along a random but isotropic direction. The 
frequency distribution of the new photon follows that of the 
source, either a Gaussian distribution with a Doppler core, 
or a continuum. Once the photon enters the gas medium, the 
length of free path is determined by calculating the optical 
depth variable r traveled during the free flight, following the 
distribution function e~ T . It is then straight forward to con- 
vert from r to distance. The location of the next scattering 
is thus determined. If the place is outside the HI cloud in 
a halo model, or if the traveled optical depth is larger than 
Gunn-Peterson optical depth in an expanding IGM model, 
the photon is labeled escaped. 

At the site of scattering, the velocity of the HI atom 
is generated by two steps. First, the velocity components 
v x and v y (normalized to Doppler velocity Vd, and z is the 
propagation direction of the photon) are generated follow- 

2 

ing a Maxwellian distribution e _ " x . Second, the velocity v z 

is generated following the distribution: f(v z ) oc ( x _^yi +a '2 
which is the joint requirement of Gaussian distribution and 
Lorentz profile for the rest frame cross section of resonant 
scattering. The direction of the resonantly scattered pho- 
tons is assumed to be isotropic, but can be easily adapted 
to other types of angular dependence. Once the direction 
is generated, frequency of the outgoing photon can be cal- 
culated. Using the notation of x, x to represent the labo- 
ratory frequency of the incoming and outgoing photon, we 
have x' — x — vcosri + vcosricos8 + vsiri6sinr]cos(f> — b(l — cos9) 



(x — V z ) +Q 



where b = m fe ^ e is the recoil parameter, r/ is the angle be- 
tween the incoming and outgoing photons, 6 and (f> are the 
two spherical coordinates of the outgoing photon where the 
coordinates are chosen such that the incoming photon is in 
z direction. (We follow Field 1959's scattering geometry and 
notations.) With this new set of frequency and direction of 
the photon, we repeat the above procedures of calculating 
the next scattering and the determination on escape. Each 
photon is followed all the way along its path until it escapes. 

Since the effectiveness of generating v z determines cru- 
cially the s peed of calculation, specia l algorit hms have been 
proposed (|Zheng fc Miralda-Escudel |2002Q, ZM02 here- 
after). We basically follow ZM02's algorithm for medium 
to large x (0.6 ^ x ^ 17). For smaller x (x < 0.6), methods 
of plain rejection (not employing ZM02's algorithm) is used. 
For very large x (x > 17), our treatment for u > ito is simi- 
lar to ZM02, but for u ^ uq, we switch the roles of the two 

_ 2 

functions, using the distribution function e " z as the trans- 
formation method to generate v z , and then use 
as the comparison function to reject. 

Twenty million of photons are experimented by Monte 
Carlo simulation for each model. Simulations are performed 
only for sources with single flash of photons. In these simula- 
tions, a time stamp can be recorded for each photon at each 
step of collision, and at its escape. For sources of arbitrary 
time dependence, a new random variable is used represent- 
ing the birth time of the photon. The time stamps can be 
generated by adding a photon's birth time to the recorded 
time stamps from a single flash simulation. Furthermore, 
with, say, 10 4 trials of randomly generated birth time for 
each original photon in a single-flash simulation, we form 
a subgroup of 10 4 photons, which have the same history of 
collisions but happen at different epochs. By this way we 
greatly improve the very low usage rate when the data of 
each recorded photon is coupled with only one birth time. 
Statistics over this enlarged group of photons give better 
continuity and smaller Poisson errors on the surveyed quan- 
tities, but do not add new physics. 



APPENDIX C: GUNN-PETERSON OPTICAL 
DEPTH 

The free path of a Lya photon in a Hubble streaming IGM 
can be derived by using GP optical depth for frequency x 
measured at source redshift z. which is 



t(x,z) = tgp(z) 



(j>(x)dx 



(CI) 



where the Gunn-Peterson optical depth parameter is given 
by equation (A3). 

In MC simulations, the Hubble streaming always makes 
frequency x lower and GP optical depth smaller. Suppose the 
photon frequency is xa immediately after the last scattering 
and xb before the next scattering, the change of GP optical 
depth St during the free flight is always negative and equals 
to -1 on average. The particular value of each St is deter- 
mined by generating a random number St ^ following the 
distribution function e Sr . With this St we can derive the 
new frequency xb by solving the equation 



St = tgp(z) 



4>(x)dx 



(C2) 
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By compiling a large data table of t(x,z) and using linear 
interpolation, xb can be solved fast. 

The free length I before the next scattering can be 
found by considering that the frequency change is caused 
by Doppler effect of Hubble motion, 



xb — xa = — 



H(z)l v 



c Ai/ D 

where Hubble constant at redshift z is 
H(z) = tf 0X /fi M (l + z) 3 + Q R (1 + zY + 

in which all Q,m, CIr, refer to values at present. 



(C3) 



(C4) 
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